// -*- C++ -*------------------------------------------------------------------
//
//  Class:              MegaAli
//
//  Base Class(es):     -
//
//  Derived Class(es):  -
//
//  Containing:         -
//
//  Authors:            Silvio Tosatto
//
//  Project Name:       BiocompUP Align
//
//  Date:               04/02
//
//  Reviewed By:        -
//
//  Description:
//    This class implements the second step in assembling a maximal subset
//    of local alignments from MegaBlast output, as required to match a genomic
//    DNA sequence against its transcript(s).
//
// ---------------------------------------------------------------------------

// Includes:
#include <MegaAli.h>
#include <IoTools.h>


// Global constants, typedefs, etc. (to avoid):

using namespace Biopool;

static void sAddLine()
{ cout << "-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-\n"; }

unsigned long int pos = 0;

inline void sResetPos()
{  pos = 0; }


inline void sDoFormat()
{
  pos++;
  if (pos % 60 == 0)
    cout << "\n";
}


inline void sWriteDot()
{
    cout << ".";
    sDoFormat();
}

inline void sWrite(char c)
{
    cout << c;
    sDoFormat();
}


// CONSTRUCTORS/DESTRUCTOR:
MegaAli::MegaAli() :  matrix(), length(), start1(), start2(), end1(), end2(), 
type1("!"), type2("!"), maxGap(10000), maxIndel(100), threshold(0.9), 
verbose(0) { }

MegaAli::MegaAli(const MegaAli& orig) 
{ 
  PRINT_NAME;
  this->copy(orig);
} 

MegaAli::~MegaAli() 
{ 
  for (unsigned int i = 0; i < matrix.size(); i++)
      matrix[i].clear();

  length.clear();
  start1.clear();
  start2.clear();
  end1.clear();
  end2.clear();

}

// PREDICATES:

void
MegaAli::writeMatrix()
{
  if (verbose < 10)
      return;

  for (unsigned long int i = 0; i < matrix.size(); i++)
  {
      for (unsigned long int j = 0; j < matrix[i].size(); j++)
	  cout << setw(3) << matrix[i][j];
      cout << "\n";
  }
}


void
MegaAli::writeStartingPoints()
{
  if (verbose <= 0)
      return;

  sAddLine();
  cout << "Possible starting points: \n";
  
  for (unsigned long int i = 0; i < start1.size(); i++)
      cout << start1[i] << "\t" << start2[i] << "\t\t" << end1[i] 
	   << "\t" << end2[i] << "\n";
}

// MODIFIERS:

void
MegaAli::copy(const MegaAli& orig)
{
  for (unsigned int i = 0; i < matrix.size(); i++)
      matrix[i].clear();

  length.clear();
  start1.clear();
  start2.clear();
  end1.clear();
  end2.clear();

  for (unsigned int i = 0; i < orig.matrix.size(); i++)
  {
      vector<long int> tmp;
      tmp.reserve(matrix[i].size());
      for (unsigned int j = 0; j < orig.matrix[i].size(); j++)
	  tmp.push_back(orig.matrix[i][j]);
      matrix.push_back(tmp);
  }

  for (unsigned int i = 0; i < orig.length.size(); i++)
      length.push_back(orig.length[i]); 
  for (unsigned int i = 0; i < orig.start1.size(); i++)
      start1.push_back(orig.start1[i]);
  for (unsigned int i = 0; i < orig.start1.size(); i++)
      start2.push_back(orig.start2[i]);
  for (unsigned int i = 0; i < orig.start1.size(); i++)
      end1.push_back(orig.end1[i]);
  for (unsigned int i = 0; i < orig.start1.size(); i++)
      end2.push_back(orig.end2[i]);

  type1 = orig.type1;
  type2 = orig.type2;
  maxGap = orig.maxGap;
  maxIndel = orig.maxIndel;
  threshold = orig.threshold;
  verbose = orig.verbose;
}

void
MegaAli::readMegaBlastData(istream& input)
{
  while (input)
    {
      string tmp1, tmp2, tmp3, tmp4, tmp5;
      long int Tlength = 0;
      long int Tstart1 = 0;
      long int Tstart2 = 0;
      long int Tend1 = 0;
      long int Tend2 = 0;
      double Teval = 0.0;
      double Tbit = 0.0;
      
      input >> tmp1 >> tmp2 >> tmp3 >> Tlength >> tmp4 >> tmp5 >> Tstart1 
	     >> Tend1 >> Tstart2 >> Tend2 >> Teval >> Tbit;

      if (!input)
	break;

      if (type1 == "!")
	type1 = tmp1;
      else
	if (type1 != tmp1)
	  ERROR("Program does not support multiple queries (yet).", exception);

      if (type2 == "!")
	type2 = tmp2;
      else
	if (type2 != tmp2)
	  ERROR("Program does not support multiple subjects (yet).", 
		exception);

      length.push_back(Tlength);
      start1.push_back(Tstart1);
      end1.push_back(Tend1);
      start2.push_back(Tstart2);
      end2.push_back(Tend2);
    }

  // find maximum values for matrix generation

  long int max1 = 0;
  for (unsigned long int i = 0; i < length.size(); i++)
    if (length[i]+start1[i] > max1)
      max1 = length[i]+start1[i];

  long int max2 = 0;
  for (unsigned long int i = 0; i < length.size(); i++)
    if (length[i]+start2[i] > max2)
      max2 = length[i]+start2[i];


  for (long int i = 0; i < max1; i++)
    {
      vector<long int> tmp;
      for (long int j = 0; j < max2; j++)
	tmp.push_back(0);
      matrix.push_back(tmp);
    }

  for (unsigned long int k = 0; k < length.size(); k++)
    for (long int i = 0; i < length[k]; i++)
      matrix[start1[k]+i][start2[k]+i] = 1;

}


void 
MegaAli::writeFragmentAli(Alignment& ali)
{
    if (verbose < 5)
	return;

   cout << "Fragment alignments:\n";
    sAddLine();
    
    for (unsigned long int i = 0; i < length.size(); i++)
    {
	
	// write fragment alignment
	
	cout << ">" << type1 << "\n";
	sResetPos();
	
	for (long int j = 0; j < start1[i]-1; j++) // unaligned
	    sWriteDot();
	for (long int j = start1[i]-1; j < end1[i]; j++)
	    sWrite(ali.getTarget()[j]);
	for (unsigned long int j = end1[i]; j < ali.getTarget().size(); j++) 
	    sWriteDot();	    // unaligned
	cout << "\n";
	
	cout << ">" << type2 << "\n";
	sResetPos();
	
	for (long int j = 0; j < start2[i]-1; j++) // unaligned
	    sWriteDot();
	for (long int j = start2[i]-1; j < end2[i]; j++)
	    sWrite(ali.getTemplate()[j]);
	for (unsigned long int j = end2[i]; j < ali.getTemplate().size(); j++)
	    sWriteDot();	    // unaligned
	cout << "\n";
	
	sAddLine();
    }
}


void 
MegaAli::concatAlignments(Alignment& ali)
{
  // init checks:

  if ( (matrix.size() == 0) || (matrix[0].size() == 0) || (start1.size() == 0)
       || (start2.size() == 0) || (ali.size() == 0) ) 
    ERROR("Invalid parameters in function call.", exception);

  cout << "Subset alignments:\n";
  sAddLine();

  vector<fragAli> fragData;
  long int bestVal = 0;

  pRecursiveSolution(fragData, bestVal);
 
  pWriteOutputData(fragData, ali, bestVal);
  writeMatrix();
  sAddLine();

  if (verbose > 1)
    {
      cout << "Control alignment:\n";
      ali.saveFasta(cout);
      sAddLine();
    }

}


void 
MegaAli::assignMaxAlignments(Alignment& ali)
{
  // init checks:

  if ( (start1.size() == 0) || (start2.size() == 0) || (ali.size() == 0) ) 
    ERROR("Invalid parameters in function call.", exception);

  cout << "Maximum subset alignments:\n";
  sAddLine();

  vector<fragAli> fragData;
  long int bestVal = 0;

  pDynamicSolution(fragData, bestVal);

  pWriteOutputData(fragData, ali, bestVal);
  sAddLine();

  if (verbose > 1)
    {
      cout << "Control alignment:\n";
      ali.saveFasta(cout);
      sAddLine();
    }
}




// OPERATORS:

// HELPERS: 

inline void
MegaAli::pSwapData(int l, int r)
{
  long int tmp = 0;

  tmp = length[l];
  length[l] = length[r];
  length[r] = tmp;

  tmp = start1[l];
  start1[l] = start1[r];
  start1[r] = tmp;

  tmp = start2[l];
  start2[l] = start2[r];
  start2[r] = tmp;

  tmp = end1[l];
  end1[l] = end1[r];
  end1[r] = tmp;

  tmp = end2[l];
  end2[l] = end2[r];
  end2[r] = tmp;
}


inline void 
MegaAli::pSortAscending(int l, int r)
{
    int i, j;
    long int v;

    if (r > l)
    {
	v = start1[r];
	i = l - 1;
	j = r;
	
	for (;;)
	{
	    while (start1[++i] < v) ;
	    while (start1[--j] > v) ;
	    if (i >= j)
		break;
	    pSwapData(i, j);
	}
	pSwapData(i, r);
	pSortAscending(l, i-1);
	pSortAscending(i+1, r);
    }
}



void 
MegaAli::pDynamicSolution(vector<fragAli>& fragData, 
long int& bestVal)
{
  // 1. sort fragments by (ascending) start1 

  pSortAscending(0, start1.size()-1);

}


inline void 
MegaAli::pRecursiveSolution(vector<fragAli>& fragData, 
long int& bestVal)
{
  // do forward recursion (i.e. initialisation)
  pForwardRecursion(fragData, bestVal);

  // do backward recursion (i.e. evaluation)
  pBackwardRecursion(fragData, bestVal);
}


inline void 
MegaAli::pForwardRecursion(vector<fragAli>& fragData, long int& bestVal)
{
  for (long int i = matrix.size()-2; i >= 0; i--)
    for (long int j =  matrix[0].size()-2; j >= 0; j--)
      if (matrix[i][j] == 0)
	continue;   // nothing to do
      else
	{
	  if (matrix[i+1][j+1] != 0)
	    matrix[i][j] = matrix[i+1][j+1] + 1; // no gap
	  else                                   // search for gapped extension
	  {
	    for (unsigned long int off = 0; off < maxIndel; off++)
	      for (unsigned long int k = 1; k < maxGap; k++)
		if ( (i+k+off < matrix.size()) 
		     && (j+off < matrix[0].size()) 
		     && (matrix[i+k+off][j]) != 0)
		{
		    if (matrix[i+k+off][j+off] > matrix[i][j])
			matrix[i][j] = matrix[i+k+off][j+off] + 1;
		}
		else if ( (i+off < matrix.size())
			  && (k+j+off < matrix[0].size()) 
			  && (matrix[i+off][j+k+off] != 0))
		    if (matrix[i+off][j+k+off] > matrix[i][j])
		    {
			matrix[i][j] = matrix[i+off][j+k+off] + 1;
		    }
	  }
	}
}


inline void 
MegaAli::pBackwardRecursion(vector<fragAli>& fragData, long int& bestVal)
{
  bestVal = 0;
  fragData.clear();

  for (unsigned long int count = 0; count < start1.size(); count++)
    {
      // do concat from start point
      
      fragAli frag;
            
      unsigned long int currI = start1[count];
      unsigned long int currJ = start2[count];
      
      frag.score = abs(matrix[currI][currJ]);

      if (matrix[currI][currJ] < 0)
	continue;

      bool valid = true;
      
      while ( (currI < matrix.size()) && (currJ < matrix[0].size()) 
	      && (valid == true) )
	{
	  aPair aP(currI, currJ);
	  frag.data.push_back(aP);

	  if (matrix[currI][currJ] > 0)
	    matrix[currI][currJ] = - matrix[currI][currJ];
	  
	  // find next pos
	  
	  if ((currI+1 < matrix.size()) && (currJ+1 < matrix[0].size()) )
	    {
	      if ((matrix[currI+1][currJ+1] != 0))    // check range
		{
		  currI++;
		  currJ++;
		}
	      else 
		{
		  unsigned long int oldI = currI;
		  unsigned long int oldJ = currJ;

		  long int maxS = 0;
		  
		  // check for gaps
		  
		  for (unsigned long int off = 0; off < maxIndel; off++)
		    for (unsigned long int k = 1; k < maxGap; k++)
		      if ( (currI+off+k < matrix.size()) 
			   && (currJ+off < matrix[0].size()) 
			   && (matrix[currI+off+k][currJ+off]) != 0)
			{
			  if (abs(matrix[currI+off+k][currJ+off]) > maxS)
			    {
			      currI = oldI+off+k;
			      maxS = abs(matrix[currI+off+k][currJ+off]); 
			    }
			}
		      else if ( (currI+off < matrix.size()) 
				&& (off+k+currJ < matrix[0].size()) 
				&& (matrix[off+currI][off+currJ+k] != 0))
			if (abs(matrix[currI+off][currJ+off+k]) > maxS)
			{
			  currJ = oldJ+off+k;
			  maxS = abs(matrix[currI+off][currJ+off+k]);

			}
		  
		  if ( (oldI == currI) && (oldJ == currJ) )
		    valid = false;
		}
	      
	    }
	  else
	    {
	      valid = false;
	    }
	  
	}
      
      fragData.push_back(frag);
      bestVal = 1;
    }

}


inline void MegaAli::pWriteOutputData(vector<fragAli>& fragData, 
Alignment& ali, double bestVal)
{
  for (unsigned int count = 0; count < fragData.size(); count++)
    if (fragData[count].score > bestVal)
      bestVal = fragData[count].score;
      
  double minValue = threshold * bestVal;

  for (unsigned int count = 0; count < fragData.size(); count++)
    {
      if (fragData[count].score < minValue)
	{
	  if (verbose)
	    cout << "Skipping short alignment.\n";
	  continue;
	}

      cout << ">>> ali# " << count << endl;
      cout << ">res1\t" << fragData[count].score << "\n";
      sResetPos();
            
      if (fragData[count].data.size() > 0)                // leading gap
	for (unsigned int i = 0; i < fragData[count].data[0].x-1; i++)
	      sWriteDot();
      for (unsigned int i = 0; i < fragData[count].data.size()-1; i++)
	{
	  sWrite(ali.getTarget()[fragData[count].data[i].x-1]);

	  for (unsigned int c = fragData[count].data[i].x-1; 
	         c < fragData[count].data[i+1].x-2; c++)
	      sWriteDot();
	}
      sWrite(ali.getTarget()[fragData[count].data[fragData[
	  count].data.size()-1].x-1]);

      if (fragData[count].data.size() > 0)                // trailing gap
	for (unsigned int i = fragData[count].data[
	    fragData[count].data.size()-1].x; i < ali.getTarget().size(); i++)
	      sWriteDot();

     cout << "\n";
     cout << ">res2\n";
     sResetPos();

      if (fragData[count].data.size() > 0)                // leading gap
	for (unsigned int i = 0; i < fragData[count].data[0].y-1; i++)
	  sWriteDot();
      for (unsigned int i = 0; i < fragData[count].data.size()-1; i++)
	{
	  sWrite(ali.getTemplate()[fragData[count].data[i].y-1]);

	  for (unsigned int c = fragData[count].data[i].y-1; 
	         c < fragData[count].data[i+1].y-2; c++)
	      sWriteDot();
	}
      sWrite(ali.getTemplate()[fragData[count].data[fragData[
	  count].data.size()-1].y-1]);

      if (fragData[count].data.size() > 0)                // trailing gap
	for (unsigned int i = fragData[count].data[fragData[
	    count].data.size()-1].y; i < ali.getTemplate().size(); i++)
	      sWriteDot();
      cout << "\n";

    }
}

